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ABSTRACT 

Observational evidence suggests a link between long duration gamma ray bursts (LGRBs) and Type Ic super- 
novae. Here, we propose a potential mechanism for Type Ic supemovae in LGRB progenitors powered solely 
by accretion energy. We present spherically-symmetric hydrodynamic simulations of the long-term accretion 
of a rotating gamma-ray burst progenitor star, a "collapsar," onto the central compact object, which we take 
to be a black hole. The simulations were carried out with the adaptive mesh refinement code FLASH in one 
spatial dimension and with rotation, an explicit shear viscosity, and convection in the mixing length theory 
approximation. Once the accretion flow becomes rotationally supported outside of the black hole, an accretion 
shock forms and traverses the stellar envelope. Energy is carried from the central geometrically thick accretion 
disk to the stellar envelope by convection. Energy losses through neutrino emission and nuclear photodisinte- 
gration are calculated but do not seem important following the rapid early drop of the accretion rate following 
circularization. We find that the shock velocity, energy, and unbound mass are sensitive to convective efficiency, 
effective viscosity, and initial stellar angular momentum. Our simulations show that given the appropriate com- 
binations of stellar and physical parameters, explosions with energies ~5x 10 50 ergs, velocities ~ 3000 kms -1 , 
and unbound material masses > 6M Q are possible in a rapidly rotating 16M Q main sequence progenitor star. 
Further work is needed to constrain the values of these parameters, to identify the likely outcomes in more 
plausible and massive LRGB progenitors, and to explore nucleosynthetic implications. 

Subject headings: accretion, accretion disks — black hole physics — gamma rays: bursts — stars: winds, 
outflows — supernovae: general 



1. INTRODUCTION 

A clear observational link has been established between 
long-duration gamma-ray bursts (LGR Bs) and Type Ic super- 
novae dGalama e t al. 1998 ] [20001 |Reichart|1999l|Bloom et al. 
20021 fDella Valle et al.||2003[ |2006[ |Garnavich et al.||2003' 



Hjorth et al.|200 3 i Kawabata et al.|2003HStanek et al.|2003 



Matheson et al. 2003; Malesani et al. 2004; Campana et al. 
2006||Mirabaletal.|2 006 ; Modjaz et al. 2006; |Pian et al.|2006[ 
Chornock et al.|2010| |Cobb et al.|2010||Starling et al.|201ip ' 
However, only a small percentage of Type Ic supernovae ex- 



hibit the late-time radio signatures of LGRBs ( Podsiadlowski 
eTaE] [20041 |Soderberg et~aLl|2"006ll. LGRBs are believed 



to be manifestations of rotationally-powered ultrarelativistic 
outflows developing in the wake of the formation of black 
holes or neutron stars in rotating progenitor. However, the 
exact mechanism for the production of LGRBs and their as- 



sociated supernovae remains a subje ct of debate (Woosley & 
|Bloom 2006 ; Hj orth & Bloom|2011| and references therein). 
At present, it is not clear whether the processes that give rise 
to LGRBs also drive a stellar explosions, or whether the ex- 
plosions are driven independently, perhaps by the standard, 
neutrino-mediated mechanism. 

In the standard supernova mechanism, an outward-moving 
shock forms after core-bounce. This shock stalls, but may 
be reinvigorated by heating by neutrinos e mitted during the 
neutro nization near the proto-neutron star ( jBethe & Wilson| 
1985), and in principle drive the star to explosion in the so- 
called "delayed neutrino mechanism." Some simulations of 
this process in at least two spatial di mensions seem to produce 
successful explosions (see, e.g., Buras et al.|2 006b; Sc heck et] 
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al 1200 6, Mezzacappa et al.|2007| [M urphy & Bu rrows|2008 
Marek & Janka 2009; Nordhaus et al.||2010b, although the 



iviaien oc jaiiK.aMz.uu y, rvuiunaus e i ai.|iz.uiu| i, aitnuugii tne 
success of two-dimensional and possibly three-dimensional 
simulations may be depend ent upon the progen i tor mass and 
the treatment of neutrinos ( |Buras et al.| 2006a; Nord haus et 
al. 2010 1. Supernovae associated with LGRBs seem to be 



more energetic than the typical Type Ic supernovae ( Iwamoto 



et al. 1998 Woosley & MacFadyen 1999 Mazzali et al. 2003 



2006b, with large kinetic energies reaching ~ 10 5 ' r ergs. Even 



if the neutrino mechanism can unbind the star, it still seems 
unclear whether it can deliver the energies found in super- 
novae associated with LGRBs. An alternative or augmenta- 
tive explosion mechanism may be required to explain the su- 
pernovae associated with LGRBs. Alternatives to the neutrino 
mechanism call on the extraction of the rotational energy of 
the central compact object — a neutron star or a black hole — 
or on tapping the gravitational energy of the material accret- 
ing toward the compact object. It remains to be determined 
which, if any, of the alternative pathways can deliver the large 
energies, and what are the resulting compact remnant masses. 

If the post-bounce neutrino-mediated energy transfer is too 
weak to unbind all of the infalling stellar strata, some mate- 
rial may continue to fall onto the proto-neutron star and possi- 
bly take it to collapse further into a black hole (e.g., Burrows 



1986, MacFadyen et al. 2001; Heger et al. 2003; Zhang et 
aU2008||Selriguchi & SMbatafioill [O'Connor & Ott|2010| ). 
This is especially relevant for rapidly rotating progenitors, 
as the progenitors with rapidly rotating cores may produce 
lower neutrino luminosities, decreasing the effectiveness of 
the neutrino-powered explosion m echanism ( |Fujimoto et al.| 
|2006l[Lee"& Ramir ez^Ruiz|2006| l. 

The infall or fallback of the stellar envelope should continue 
past the initial emergence of the event horizon, but then the 
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structure of the accretion flow becomes sensitive to its angu- 
lar momentum content. Given sufficient angular momentum, 
the flow becomes rotationally supported. Such a "collapsar" 
configuration has been proposed to naturally lead to the ultra- 
relativistic outflow in LRGBs ( Woosley 1993 ), as gamma rays 
can be produced in an ultrarelativistic jet launching from the 
magnetosphere of the black hole that forms in the aftermath 
of the collapse of the rotating progenitor. The jet is powered 
by a continuous infall and disklike accretion of the progenitor 
star's interior. 

It has long been hypothesized that a "wind" outflowing 
from a collapsar accretion disk could unbind the stellar en- 
velope and synthesize su fficient 56 Ni to produce an optically 
bri ght supernova (e.g., MacFadyen & Woosley|199 9 1 Pruet et 



(Milosavlj evic et aT||2012 i suggest that followi ng shock for 



aLl[2003| [2004} |Kohri et al.||2005| l. The dynamics of the en 
ergy flow in such a system has yet to be elucidated. In the 
present study, we utilize one dimensional hydrodynamic sim- 
ulations with rotation (1.5D) to test the hypothesis that ac- 
cretion power can drive an explosion of the star. We do not 
simulate the core bounce, and simply posit that any prompt 
and neutrino-reinvigorated shock has failed and that the stel- 
lar atmosphere has not acquired outward motion and is free to 
ac cret e toward the black hole . 

In |Lindner et al.| ( 2010| l, we simulated the post-core- 
collapse hydrodynarmcal evolution of a rapidly r otating 
14M Q Wolf-Rayet stellar model 16TI of (jWoosley & Heger 



2006) that has been proposed as a LGRB progenitor. The 



rate at which the infalling stellar envelope was being accreted 
onto the black hole evolved through two distinct phases dur- 
ing the first ~ 200 s following the initial collapse of the stellar 
core. First, the low specific angular momentum material of 
the inner layers of the star accreted quasispherically through 
the inner boundary, and is presumed to have accreted onto 
the black hole. Then, the material that had sufficient angular 
momentum to become rotationally supported on the computa- 
tional grid formed a thick accretion torus. Simultaneously, an 
accretion shock appeared at the innermost radii and traversed 
the star. Most of the stellar envelope traversed by the shock 
was in radial hydrostatic equilibrium and convective; convec- 
tion transported the energy dissipated at the smallest simu- 
lated radii toward the expanding shock. The central accretion 
rate was nearly time-independent prior to rotating torus and 
shock formation, and dropped sharply afterwards. The abrupt 
drop of the accretion rate closely resembled the prompt 7- 
ray and the early X-ray LGRB light curves mea sured with 
the NASA Swift satellit e (Tagliafe rn et al.||2005| INousek et] 
|aL|2006[ IQ'Brien et al.||2006| >, adding weight to the hypofh- 
esis that the light curves are responding to an e volution of 
the central accretion rate ( |Kumar et al . 2008a b). Because 
the innermost simulated radius was 500km, much larger than 
the innermost stable circular orbit around the central black 
hole (5 -50km), the accreted-mass-to-energy conversion effi- 
ciency was low and the shock acquired relatively low veloci- 
ties, ~ lOOOkms" 1 , while in the interior of the star. The star 
did not explode, but only lost mass to the thermally driven 
wind that set in after the shock had traversed the star. 

In collapsars, a substantially larger accretio n energy is dis- 
sipated at the radii left out from the [Lindner et al.| ( |2010| ) 
simulations, closer to the black hole, but only a fraction of 
this energy couples to the stellar envelope. The rest may be 
lost to the emission of neutrinos and to the photodisintegra- 
tion of hydrostatic elements into free nucleons as well as to 
advection into the black hole. Crude analytical considerations 



mation and the rapid accretion rate drop seen in Lindner et al. 
(2010), neutrino losses are relatively small. Then, the amount 



of energy transferred onto the envelope is determined by the 
competition of the inward advective and the outward convec- 
tive energy transport. The advection arises from the inward 
drift of the fluid in response to magnetohydrodynamic (MHD) 
stresses; the convection arises from entropy gradients arising 
from the dissipation of MHD turbulence. If convective trans- 
port is efficient, the amount of energy transferred from near 
the black hole to the shocked envelope can be sufficient to 
drive a fast shock with velocity ^ 1000km s" 1 a nd unbind 
the star. The model of Milosavljevic et al. (2012 1 suggests 



that the parameters determining the viability and energy of 
such accretion-powered supernovae are the viscous stress-to- 
pressure ratio a and the convective mixing length A conv . The 
model could not, of course, capture the consequences of the 
interplay of pressure and rotation at the critical radii where the 
two sources of radial support against gravity are comparable. 

In this work we show the results of a series of rotating 
one-dimensional simulations of the immediate aftermath of 
the collapse of a rapidly rotating LGRB progenitor star's 
core. While one-dimensional, our simulations include rota- 
tion in a spherically-averaged sense and implement a modi- 
fied a-viscosity prescription. One customarily refers to such 
simulations as "1.5 dimensional." They also take into ac- 
count optically thin cooling by neutrino emission, cooling and 
heating by nuclear processes, and energy and compositional 
transport by convection in the mixing length theory approx- 
imation. This work is complementary to our rotating two- 
dimensional sim ulations (2.5D) of collapsar accretion ( jLind- 



ner et al. 2010), in which we simulated only relatively large 



radii and did not incorporate nuclear and neutrino physics. 
Here, we sacrifice in spatial dimensionality to make it pos- 
sible to track rudimentary nuclear compositional transforma- 
tion and simulate smaller radii (r > 25 km) over similarly ex- 
tended time periods (~ 40- 100s). In the presence of cooling 
by neutrino emission the rota ting central torus may be geo- 
metrically thin (e.g., MacFadyen & Woosley||1999| [Popham 



1999 1 IKohri et al.||2005t |Chen & Beloborodov PoTT 
guchi & Shibata 201 l||Taylor et al.|2010|). Therefore, we 



m 

'7; 



include corrections to approximate the effects of such flow. 
The principal source of model uncertainty is the efficiency of 
convection which in the mixing length approximation can be 
parameterized with an effective value of the mixing length. 
To our best knowledge, there has not been a systematic first 
principles study of convective efficiencies in the rapidly con- 
vecting regime. Thus the mixing length A con v, along with the 
viscous shear stress-to-pressure ratio a, are the parameter de- 
pendences that we explore. 

A magnetic outflow driven by a proto-neutron star may 
carry an energy similar to that of a supernova (e.g., 
Bisnovatyi-Kogan 1971; Wheeler et al. 2000; Thompson et 
al.|2004HBucciantini et al.|20071|Burrows et aL|20071|Dessart 



et al. 



2008 ). However, the outflow may be too axially colli- 



mated to produce a standard, quasi-spherical explosion ( Buc- 



ciantini et al.|2008 , 2009 ). Here, we assume that any explo- 



sion mechanism preceding the collapse into a black hole has 
failed. Clearly, our one-dimensional model cannot capture the 
effects of the formation of a magnetized jet, after an accretion 
disk has formed. Although this is an integral component to 
the collapsar model for LGRBs, we omit any treatment of the 
jet in the present work. 
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This work is organized as follows. In Section [2] we dis- 
cuss our numerical algorithm. In Section [51 we present the 
results our simulations. In Section [4j we identify the param- 
eters critical to our model, and discuss their implications for 
real accretion powered supernovae. Finally, in Section [5] we 
summarize our conclusions. 

2. NUMERICAL ALGORITHM 

The simulations were carried out with the piecewise- 
parabolic method (PP M) solver in the adaptive-mesh- 
refinement code FLASH (Fry xell et al.|2000[ ), version 3.2, in 
one spatial dimension. Although the rotating stellar collapse 
is inherently three-dimensional, we have chosen to approxi- 
mate the key multidimensional effects, including angular mo- 
mentum transport and convective energy and compositional 
transpo rt, w ith a spherically-averaged transport scheme. In 
Section |2~T| we describe our imp lementation of angular mo- 
mentum transport. In Section 2.2 we describe our calculation 



of the self-gravity of the fluid. In Section 2.3 we describe our 



modeling of the transition toward nuclear statistical equilib- 
rium (NSE) in the hot inner accretion flow. In Section[Z4] we 
discuss cooling by neutrino emission. In Section [23] we de- 
scribe our treatment of conve ctive energy transport and com- 
positional mixing. In Section 2.6 we describe the corrections 
that we apply in situations where, in the presence of cooling, 
the acc retio n flow is expected to be geometrically thin. In 
Section |2~7| w e describe our initial and boundary conditions. 
In Section |2T8l we sh ow the results from tests of the code. Fi- 
nally, in Section[Z9j we briefly review the various limitations 
our method. 

2.1. Angular Momentum 

To include rotation and angular momentum transport in our 
one dimensional model, we track the specific angular momen- 
tum £ = rv^, where is the azimuthal velocity, which we in- 
terpret as the mass-weighted spherical average of an underly- 
ing polar-angle-dependent angular momentum £(r,8). If, e.g., 
spherical shells rotate rigidly, £(r, 9) oc sin 2 #, and the fluid 
density is spherically symmetric, then the one-dimensional 
specific angular momentum is two-thirds of the midplane 
value, i = |£mid. The azimuthal Navier-Stokes equation, com- 
bined with the equation of mass continuity, then implies the 
one-dimen sional angular m omentum transport equation (see, 
e.g., [Thompso n et al.|2005| ) 



d(p£) + 1 <9(r 2 v r/ c 



dt 



Or 



1 9 I 3 \ 



= 0, 



where v is a shear viscosity and 



d 



dr 



(1) 



(2) 



is the r— <p component of the shear tensor. The energy dissi- 
pated through shear viscosity was accou nted for by includi ng 
the specific heating rate (see, e.g.,|Landau & Lifshitz|1959]l 



Q 



vise 2 

= va rrh 



(3) 



where gvisc denotes the volumetric viscous heating rate. The 
dimensional reduction in equation ([T]) is inaccurate in re- 
gions where the disk is geometrically thin. There, the mass- 
weighted spherical average closely approximates the mid- 
plane value, £ ~ l^. We ignore this effect, but we do incor- 



porate thermodynami c co rrections addressing the transition to 
a thin disk in Section [2~6l 

Our tr eatment of shea r viscosity is similar to our method- 
ology in |Lindner et al.| ( |2010| l, and for completeness we re- 
produce our methodology here. Since we do not simulate the 
magnetic field of the fluid, we utilize a local definition of the 
shear viscosity to emulate the magnetic stress arising from 
the nonline ar development of the m agnetorotational instabil- 
ity (MRI; Balbus & Hawley 1998 and references therein). It 
should be kept in mind, however, that the effects of MRI are in 
some respects very different from those of the viscous stress. 
For example, the thick disk surrounding our collapsar black 
hole is convective; in unmagnetized accretion flows convec- 
tion transports angular momentum inward, toward the center 



of rotation (Ryu & Goodman 1992 Stone & Balbus 1996 



Igumenshchev et al. 2000 ), whereas in magnetized flows, con 



vection can also transport angular momentum outward (Bal- 



bus & Haw ley 2002 ; Igumenshchev 2002) |Igumenshchev et 

Although we include 



al.| |2003| [Christodoulou et al.||2003| > 
treatment for co nvec tive energy flux and compositional mix 
ing (see Section 2.5 1, we do not include angular momentum 
transport by convection. 

Our definition of the local viscous stress emulating the 
MRI must be valid under rotationally supported, pressure 
suppo rted, and freely fa lli ng conditions, and we proceed as 
_ d20T0k. 
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Thompson et al. (20051 suggest 
that since the wavenumber of the fastest growing MRI mode, 
which is given by the dispersion relation vpjt ~ fl where va 
is the Alfven velocity and Q = v^/r is the angular velocity, 
should in the saturated quasi-state state be about the gas pres- 
sure scale height, k oc H~ l , the Maxwell pv\ and viscous vp^l 
stresses (up to factors in |c/lnil/oflnr| that we neglect) can be 
equated if the viscosity is given by 



^mri = aH fl, 



(4) 



where a is a dimensionless parameter. If the pressure scale 
height is defined locally, 

~\ (5) 



H = iVlnPl 



the viscosity defined in equation (RJ) suffers from dive rgences 
at pressure e xtrema. To alleviate this problem, as in |Lindner| 
|et al.| (|2010[), we define a second viscosity according to the 
Shakura & Sunyaev ( 1973 1 prescription 



P 



(6) 



Shakura-Sunyaev viscosity overestimates the magnetic stress 
in stratified hydrostatic atmospheres. We thus set the viscosity 
in equations ([T]i and Q to equal the harmonic mean of the 
above two viscosities 



2 z^mri vss 
^mri + vss 



(7) 



where the pressure gradient in equation d5]l is calculated by 
the finite differencing of pressure in neighboring fluid cells. 
Additionally, we have applied a Gaussian kernel smoothing 
to the radial dependence of H to help filter short-wavelength 
num erical instabilities. We describe this procedure in Section 
E3J 

In FLASH, we treat specific angular momentum as a con- 
served "mass scalar" that is being advected with the fluid, 
which makes p£ a conserved variable; the corresponding cen- 
trifugal force is then incorporated in the calculation of the 
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gravitational acceleration, as we explain in Section 2.2 below. 
Then the third parabolic term in equation ([T]) is computed ex- 
plicitly through the inclusion of the radial p^-flux -rvpa r( j, in 
the advection of I. 

Numerical stability of an explicit treatment of a parabolic 
term places a upper limit on the time step 



self gravity is calculated from 



At< 



Ar2 



(8) 



where Ar is the grid resolution. For a 0.01, the viscous 
time step in our simulations becomes significantly shorter 
than the Courant time st ep. In our test inte grations with a 
7-law equation of state (Lindner et al. 2010[ ), we find that, 
while not implying an outright instability, a choice of Af that 
saturates the limit in equation ([SJl results in weak stationary 
staggered perturbations in the fluid variables. We ignore this 
complication and allow our time step to be set by the limit 
in equation ([8]l of the cell with the smallest viscous diffusion 
time across the cell. 



2.2. Gravity 

We calculate contributions to the gravitational potential 
from a central point mass and a spherically-symmetric ex- 
tended envelope. General relativistic effects become impor- 
tant at the innermost radius, which in some simulations is as 
small as r m i n = 25 km. At radii r ~ r^, the black hole dom- 
inates the enclosed mass after about 0.5 s. Thus, we describe 
the gravity of the black hole using the approximate, pseudo- 
Newtonian gravitati onal force for a rotating black hole pro- 
posed by Artemova et al.j (|1996| ), which is a generalization 
of the Paczyhski & Wiita[( [l980[ ) pseudopotential to rotating 
black holes. However, we continue to calculate the gravity of 
the fluid in the Newtonian limit. The Artemova et al. gravita- 
tional acceleration in the equatorial plane of a rotating black 
hole is given by 



gBH(r,0 = ir/2) = - 



GM B 



- /3 (r-r H ) /3 



(9) 



where re = [1 +(1 -a^'^IGMbh/c 2 is the radius of the event 
horizon expressed in terms of the dimensionless spin param- 
eter a, and f3 = r^co /rn — 1 is a dimensionless exponent with 
risco denoting radius of the innermost stable prograde equato- 
rial circular orbit. We assume a dimensionless spin parameter 
of a = 0.9 in these calculations. Our treatment does not incor- 
porate general relativistic corrections to the viscous stress and 



momentum equations (see, e.g., Beloborodov| 1999 1. 

We adopt the form of the gravitational acceleration in equa- 
tion |9]), which was derived for the equatorial plane of the 
black nole, to represent the mass-weighted spherical average 
of the gravitational acceleration, by setting gBn(r) = gsuif, 9 = 
7r/2). This approximation is appropriate when the accret- 
ing mass is concentrated in the equatorial plane, especially 
when the innermost disk is geometrically thin, and is prob- 
ably rather inaccurate for an accretion flow that is geometri- 
cally thick down to r^co- Our simulations predict a geomet- 
rically thin disk at r < 100 km or greater radii after material 
has circularized in our simulation, so this assumption seems 
adequate. 

For each zone, the gravitational acceleration due to fluid 



AttG J 



3 / & r < 



ArA' 



ArA' 



■?(10) 



where Ar, and Ar^ are the radial widths of the grid cells. The 
net gravitational and inertial acceleration in our calculation is 
then given by 

a t ot = gBH + gself + a C ent , ( 1 1 ) 



where 



Scent — ~~ 5~ r 



(12) 



is the centrifugal acceleration. 

2.3. Nuclear Processes and the Equation of State 

To calculate the internal energy of t he fluid, we use the 
Helmh oltz equation of state (EOS) of |Timmes & Swesty 
(2000) included with the FLASH distribution, which ac- 
counts for the contributions to pressure and other thermo- 
dynamic quantities from radiation, ions, electrons, positrons, 
and Coulomb corrections. We track the abundances of 47 nu- 



clear isotopes treated in the NSE calculations of Seitenzahl 
et al. (20081 and pass the local nuclear composition to the 
EOS as input. Given density, temperature, and nuclear com- 
position, the Helmholtz EOS provides the internal energy, 
density, pressure, entropy, specific heats, adiabatic indices, 
electron chemical potential, and various derivative thermody- 
namic quantities. During the course of the thermodynamic 
update and the cooling update which is operator-split from 
the thermodynamic update, the temperature must be derived 
from the internal energy, and in the Helmholtz EOS, this is 
achieved by numerically solving for the implicit relation 



£EOs(p,r,X) = e 



(13) 



for the temperature, where e is the specific internal energy and 
X = (Xx , ...,Xn) is the vector of isotopic mass fractions X t . 

The fluid heats and cools in response to nuclear composi- 
tional transformation. We do not integrate a nuclear reaction 
network, but instead model the change of the nuclear com- 
position as a gradual convergence to nuclear statistical equi- 
librium (NSE) in the part of the flow where the convergence 
time scale t N se is comparable to or shorter than the age of 
the system. In this model, as we explain below, the nuclear 
composition responds instantaneously to a change of the tem- 
perature, implying that the dependence of the composition on 
the temperature must be taken into account, in a manner that 
conserves the combined specific internal and nuclear energy 
e + e nuc when solving the EOS for temperature. Here, e nuc is 
the specific (negative) nuclear binding energy of the fluid 



4^ A,m„ 



(14) 



while E B i is the negative nuclear binding energy of the isotope 
and Aj is the atomic mass of the isotope. 

The time scale for conv ergence t o NSE can be app roxi- 
mated via ( |Khokhlov|199T1 see, also, |Calder et al.|2007| ) 



/ 279 7 

tnse = p°' 2 exp [—^-- 40 - 5 1 S ' 



(15) 
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where T = 10 9 raK and p is the density in gem" 3 . At rele- 
vant densities, this time scale is of the order of one second 
for 7nse ~ 4 x 10 9 K. We calculate the nu clear mass frac- 
tions u sing the publicly available solver of Seitenzah l et al.| 
(2008) which solves for the NSE mass fractions X^se./ of 47 
nuclear isotopes as a function of density p, temperature T, 
and proton-to-nucleon ratio Y e = J^-Z^/A;, where Z, is the 
atomic number an isotope. At temperatures T > 3 x 10 9 K we 
model convergence to NSE via 



dXi_ 
~di 



^ r ;,NSE(p, Tnsei Y e ) —Xj 
tnse(p, T) 



(16) 



where T N se is the temperature that the fluid element would 
have given enough time to relax into NSE while keeping the 
total specific energy e + £nse and proton-to-nucleon ratio Y e 
fixed. The temperature Tnse is implicitly defined by the con- 
dition (cf. equation (13}) 

e EOs[/3, 7nse,X N se(P, Tnse,^)] +£nuc[XNSE(P, T N s E ,Y e )] 

= e(p,r,X) + e nuc (X). (17) 

This condition ensures that the sum of the internal and nu- 
clear energy densities in NSE would equal the sum of the two 
energy densities in the model. We solve equation ( [T7| for 
T N sE(p,Y e ,e,e nuc ) iteratively and then update the abundances 
by discretizing equation ( 16 1 with 



X i (t+At) = X i , NSE (p,T NSE ,Y e ) 

+[Xi(t)-X imE (p, r NSE , Y e )] exp 



At 



tnse(p, T) 



(18) 



Following the update of the nuclear mass fractions, we update 
the specific internal energy to account for heating or cooling 
due to any change in specific nuclear binding energy 



e(t + Af ) = e(t) + e nuc [X(f )] - e n uc [X(f + Af)] , 



(19) 



and finally update the temperature from equation ( 13 i 



This prescription does not affect the proton-to-nucleon ra- 
tio Y e ; that latter is a conserved mass scalar in our simulations. 
Thus, the expected partial neutronization in the mildly degen- 
erate innermost segment of the accretion flow is not calculated 
and our prescription cannot be used to accurately estimate the 
56 Ni fraction within the Fe-group elements synthesized in the 
simulation. 

2.4. Cooling 

The hot innermost accretion flow cools via neutrino emis- 
sion. At the densities observed in our simulation, the disk and 
stellar atmosphere are transparent to neutrinos. The two most 



significant neutrino-emission channels (e.g., Di Matteo et al. 
|2002| and references therein) are: 

1. Pair capture on free nucleons (the Urea process): p + 
e~ n + v and n + e + —> p+v. The cooling rate is 



Q eN = 9 x lO^pio^n^nuc ergs cm 3 s 



(20) 



where p= lO'Viogcm -3 , T = 10 n T n K, andX nuc =X p +X„ is 
the mass fraction in free nucleons. 
2. Pair annihilation (e~ + e + — > v + v). The cooling rate is 



Q e+e - = 1 .5 x KPT?! ergs cm" 3 s _1 . (21) 
All three flavors, e, p, and r, of neutrinos are included. 



We have included the above neutrino cooling rates in our 
calculations, where losses are computed via 



e(f + Af) = e(f)-— Af, 
P 



(22) 



where Q u = Q e n + Q e + e - is the total volumetric neutrino cool- 
ing rate. The update of the internal energy due to cooling is 
operator split from the update due to nuclear compositional 
change. 

2.5. Convection 

We introduce convective energy transport and composi- 
tional mixing within the framework of mixing length theory 



(e.g.,|Kuhfu6 1986). In the calculation of the convective trans- 
port fluxes, we ignore the radial variation of the mean molecu- 
lar weight as well as rotation, and the condition for instability 
is simply the Schwarzschild criterion, ds/dr < 0. Then, in 
unstable zones, the convective energy flux is 



1 



conv — ^ CppVcon V : i 



1? A c 



dT 

a7 



ds 
dr ' 



(23) 



where c P is the specific heat at constant pressure, A conv is the 
length over which convection occurs, s is specific entropy, and 
Vconv is the convective velocity. The convective velocity can 
be approximated by 



1 



r An 



dp 
df 



T ds 
P cp dr 



1/2 



(24) 



where g < is the gravitational acceleration in the local rest 
frame of the convectively unstable fluid 



dv 
dt 



1 dP 
p dr 



P 

pH> 



(25) 



and g erav = gBH+gseif is the net gravitational acceleration in the 
inertial frame, v in the second step denotes the mass-weighted 
spherical average of the fluid velocity at radius r. To pa- 
rameterize our uncertainty regarding the value of the convec- 
tive mixing length, we introduce a dimensionless parameter 
Cconv ~ 0(1) defined as 



Then, combining equations ( [23 
expression 



conv 



(26) 



261, we obtain the standard 



^conv — ^sconv^ Cp 



pfdp 

H\dT 



1/2 



T ds\ V2 
cp dr J 



(27) 



which is appropriate even when the fluid is not in hydrostatic 
equilibrium and v r ^ 0. 

In evaluating the convective energy flux at a boundary 
(face) of a computational cell, we use face-centered linear 
interpolation of the density, temperature, and pressure. The 
zone-centered values of the specific heat c P , specific entropy 
s, and thermodynamic derivatives {8P/dT) p , and (dP/dp)r 
are returned by the EOS routine, and the face-centered values 
are again computed by linear interpolation. Then, (dp/dT)p 
is calculated from 



dp 
dT 



OP 

df 



dP 

dp 



(28) 
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The convective energy flux never exceeds 

^conv < pec s , (29) 

where c s = ("fcP/p) 1 ^ 2 is the adiabatic sound speed, and 7 C is 
the adiabatic index. 

We anticipate that a local application of MET, in which the 
expression for the convective energy flux contains a pressure 
derivative in the denominator, may contain an instability. The 
instability is an artifact of modeling the intrinsically nonlocal 
convective energy transport with a local nonlinear differen- 
tial operator. To control — if not entirely prevent — undesirable 
outcomes of the instability, we filter short wavelength pertur- 
bations in the calculation of the pressure scale height H that 
enters our estimates of the viscosity and the energy flux trans- 
ported by convection by applying a Gaussian smoothing 



^smoothC^") " 



(30) 



where the summations are over all of the cells in the sim- 
ulation, and the spherically-averaged smoothing kernel is 
given by 



k i (r) = 



1 Ann 



2V2tt ra 



exp 



(r-n? 
2a 2 



-exp 



(r+nf 
2a 2 



(31) 

Here, a is a radius-dependent smoothing length that we set to 
\r. Similarly, in the evaluation of the specific entropy deriva- 
tive in equation p7), we smooth the specific entropy s via 



■^smooth \X) — 



Y,M r ) Pi s i 



(32) 



The filtering affects only the evaluation of F com and helps 
avoid breakdown of our transport scheme, but residual ar- 
tificial non-propagating waves do develop, and saturate, on 
wavelengths comparable to the smoothing length. 

The accretion shock formally presents a negative entropy 
gradient but physically does not give rise to convection. The 
upstream of the Shockwave is marginally convectively stable 
as the Shockwave traverses the progenitor's convective core, 
and becomes absolutely stable in the radiative envelope. To 
prevent spurious convection across the shock transition, we 
modify the convective flux to decline to zero linearly near the 
shock 



i(r) : 



(l-^Ashock^convW, r < r shock, 



0. 



Y f shock i 



(33) 



where r s h oc k is the radius of the accretion shock front which 
w e track during the si mulation. 

|Murphy & Meakin] ( |201 \\ argue that on physical grounds, 
in quasi- stationary "stalled" shocks in the standard core- 
collapse context, the distance from the shock r s h oc k — r is the 
appropriate convective length scale near the shock, as convec- 
tive eddies can grow to the largest size available to them. If we 
had set the convective mixing length A conv proportional to the 
distance from the shock, which is the adaptation of MET that 
Murphy & Meakin suggest, equation p3] > would have con- 
tained a quadratic factor (1 - r/r s h 0c k) 2 , instead of the linear 
factor (1 -r/r s h oc k) that we employ. The physically motivated 
modification of A conv of Murphy & Meakin, which we became 
aware of after the completion of this work, and our ad hoc ver- 
sion should give rise to similar dynamics, especially when the 
shock travels outward as in our simulations. 



Convection also gives rise to compositional mixing in the 
convective region. We model the mixi ng of nuclear s pecies 
in the diffusio n approximation (e.g., Cloutman & Eol l|1976 
IKuhfuBI 1986)1 



d(pXj) 
dt 



]_d_ 

r 2 Q r 



(f J~ mix,/) i 



where 



I 9* 

3 Or 



(34) 



(35) 



is the mass flux of species i transported by convection, while 
t'conv is the compositional diffusivity which we take to be pro- 
portional to the convective velocity multiplied by the pressure 
scale height 

^conv — £mix ^conv Aconv 5 (36) 



and £ m ; x ~ 0(1) is a dimensionless parameter. We again apply 
the flux limitation behind the shock front in the form of the 
linear factor in equation ( [33) . The compositional diffusion 
is also subject to the timestep limitation imposed in equation 
"S). It is worth noting that compositional diffusion implies a 
ux of nuclear energy given by 



AiTTlp 



(37) 



The entropy transport equation implied by our algorithm is 

■Qu + Qmo, (38) 

where d/dt = d/dt +v r d/dr and 



ds 1 d 2 jp \ n 

pi ~T ~l 9 "Ti - V ^conv.mod) — \dv 

at r L or 



Gnuc = ~P 



Ajtrip dt 



(39) 



is the rate of heating or cooling associated with nuclear com- 
positional transformation (see equations 1 14 1 and 1 16 1). 



2.6. Thin Disk Corrections 

We have thus far assumed a rotating, quasi-spherical ac- 
cretion flow. However, near the black hole, where cooling 
by neutrino emission and nuclear photodisintegration into nu- 
cleons is significant, the flow can becom e geome trically thin 
(e.g.,[M acFadyen & Woos ley 1 1 999] [ Popham e t al |1999|[Kohri] 
|et al.)|2005| |Chen & Beloborodovl|2007| >. In this case, the 
quasi-spherical treatment underestimates the density, pres- 
sure, and temperature near the midplane of thin disk, where 
the bulk of the neutrino emission takes place. We introduce a 
correction that adjusts the temperature of the flow to be closer 
to the physical, thin-disk value. 

In what follows, the quantities applying to the thin disk will 
be marked with tilde. Let H z denote the vertical half-thickness 
of the thin disk. We assume that the vertical half-thickness 
of the quasi-spherical flow is H z ~ |r. Since the thin and 
the quasi-spherical flow must contain the same column den- 
sity, H,p ~ \rp. Ignoring differences in nuclear composition 
between the thin and quasi-spherical flows, the same corre- 
spondence must apply to the total internal energies integrated 
along the vertical column, H z pe ~ \ r P^ an d thus, e ~ e while 
H Z P ~ ^rP. Vertical force balance in the thin disk requires 

P/H z = p\g z \, where g z = -sgn(z)gH z /r is the gravitational 
acceleration in the z-direction, and g = (g B H + gseif) • r is the 
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radial gravitational acceleration (see Section 2.2 1. Enforcing 
that H- < H z , we obtain 



H 7 = min 



rP\' r 
P8 



(40) 



To account for the higher density in the thin disk, we could 
pass p to the EOS. This, however, would result in a modified 
pressure P. Out of a possibly unfounded concern that a mod- 
ification of the pressure would introduce spurious dynamics 
in the spherically-averaged flow, we opted to modify our es- 
timate of the disk midplane temperature in a manner not di- 
rectly affecting the fluid pressure. This corrected temperature 
then enters the calculation of the neutrino cooling rate and the 
NSE composition, both of which are highly sensitive to the 
midplane temperature. 

We estimate the midplane temperature T from the following 
extrapolation 



dlnT 
dlnp 



Ml" 

P 



(41) 



where the partial derivative, which we denote with x, is evalu- 
ated at constant specific internal energy and can be expressed 



ainr 

dlnp 



>-(—) 
\d\n P J T 



ds 



dlnT J pc v T 



(42) 



where cy is the specific heat at constant volume. The quanti- 
ties on the right hand side of equation ( |42| , are all provided 
by the Helmholtz EOS. 

Since p/p ~ r/(2H,), the midplane temperature of the disk 
can be approximated via 



r 
2H, 



T = ET, 



(43) 



where the last equality defines the dimensionless temperature 
correction factor S. To ensure continuity near the shock tran- 
sition, we modify the correction factor to linearly approach 
unity at the shock transition by defining 



^mod " 



1+1 



F shock 



(S-l). 



(44) 



For clarity of notation, we drop the subscript in S mo d in what 
follows. 



The correction introduced in equation ( 43 1 affects both the 
temperature calculated from internal energy via the equation 
of state, as well as the NSE temperature calculated from the 
total internal and nuclear energy as described in Section |2.3 
above. Equation ( 13 1 is corrected to become 



£eos(P,S 1 f,X) = e, 
while equation ( fT7| is corrected to become 

£eos[p,S 3nse,Xnse(Pj TksE,Y e )] 

+ £nuc[XNSE(P, 7NSE,ie)] 

= e( j0 ,S- 1 f,X) + e nuc (X). 



(45) 



Note that since H > 1, the estimated midplane disk temper- 
ature is higher than the temperature calculated without this 
correction, but this allows the disk to cool faster than it would 
otherwise. The rate of cooling by neutrino emission is then 
calculated from equations p0| ) and pT) but at density p and 
temperature T. 

2.7. Initial Model and Boundary Conditions 

The initial model is the ro tating M stal « 14M Q Wolf-Rayet 
star 16TI of Woosley & Heger (2006), evolved topre-core- 
collapse from a 16M© main sequence progenitorF To pre- 
pare the model 16TI, Woosley & Heger assumed that the 
rapidly rotating progenitor, which is near breakup at its sur- 
face at r star w 4 x 10 s km, had low initial metallicity, 0.01 Z Q , 
and became a WR star shortly after central hydrogen deple- 
tion, which implied an unusually small amount of mass loss. 
For illustration, the specific angular momentum at the three- 



quarters mass radius was 



■3/4 ' 



-8x10 cm s , implying cir- 



cularization around a 5M black hole at r ~ 2500km, much 
larger than ISCO. The circularization radii of the outermost 
layers of the star are in the range 10 4 - 10 5 km. Woosley & 
Heger provide a radius-dependent angular momentum profile 
£\6Ti(r)- We introduce the dimensionless parameter ^ to scale 
the specific angular momentum £(r) of our initial model rela- 
tive to that of 16TI 



«r) = £<Wr). 



(47) 



The plots of density, temperature, angular momentum, and 
composition in Section|3|show the initial conditions. The an- 
gular momentum profile is specific to our fiducial Run 1 with 
£e = 0.5, half of the rotation rate of 16TI. 

The iron core of the model 16TI, with a mass ~ 1M , has 
mass too low to collapse directly into a black hole, but should 
instead first collapse into a neutron star. The latter could, but 
need not to, be driven to a successful explosion by the delayed 
neutrino mechanism. A black hole can form by fallback. We 
do not in any way account for the core bounce and its con- 
sequences, nor for the heating by the neutrinos emitted from 
the proto neutron star. Our central compact object is a point 
mass from the outset equipped with, as we clarify below, an 
absorbing boundary condition. 

Pseudo-logarithmic gridding is achieved by capping the 
adaptive resolution at radius r with Ar > \rir where r\ is a 
dimensionless parameter. We choose rj = 0.15 for all but Run 
2, where rj = 0.075. Beyond the outer edge of the star we 
place a cold (10 4 K), low-density, stellar- wind-like medium 
with density profile p(r) = 3 x 10~ 7 (r/r stal )" 2 g cm" 3 . 

The simulation was carried out in the spherical domain 
'"min <r < r max . We placed the inner boundary at r m ; n ~ 25 km 
and the outer boundary well outside the star at r max = 10 7 km. 
In Table [T] we summarize the main parameters of our simula- 
tions, and also present some of the key measurements, defined 
in Section [3] characterizing the outcome of each simulation. 
Each simulation was run for ~ 10 7 hydrodynamic time steps 
and required ~ 5000 CPU hours to complete. 

The boundary condition at r lmn was unidirectional "out- 
flow" that allowed free flow from larger to smaller radii (v r < 



(46) 



1 |L6pez-Camara et~al l j2009t carried out SPH simulations of neutrino- 
cooled accretion during the first 0.5 s of the collapse and [Morsony et aT] 
1 2007 1 and Nagakura et al. ( 20 1 1 a I simulated the propagation of a relativis- 
tic jet using the same model star. We discuss important caveats of using this 
model in Section|4] 
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FIG. 1 . — Test of internal energy conservation in a run identical to Run 2 except with nuclear compositional change disabled and using a Newtonian gravitational 
potential. Plotted are each of the terms in equation \50) calculated at r m j„ = 50km (left) and r m j„ = 1000km (right). Note the difference in scales on the vertical 
axis. The quantity AE represents the absolute value of the difference between the left and right hand side or equation \5Q\ . A possibly dominant source of 
apparent error is inconsistent discretization of the various fluid variables and fluxes in PPM and in the post-processing and need not reflect an inaccuracy of the 
computation. Ne utrin o losses are insignificant outside the inner ~ 1000km. Over the interval 0s < t < 70s, we find that the total error is < 1% of the largest 
term in equation {50). 



0) and disallowed flow from smaller to larger radii (v r > 0) 
by imposing a reflecting boundary condition. We imposed the 
torque-fr ee boundary condition via (see, e.g., Zimmerman et 
|al.|2005D 

m, - 

As in other Eulerian codes, the boundary conditions in 
FLASH are set by assigning values to fluid variables in rows 
of "guard" cells just outside the boundary of the simulated 
domain. Let r\n denote the leftmost cell within the simu- 
lated domain, and let rg where Q = (— I,— | ,— I ,— |) be the 
four guard cells to the left of r\/2 such that the grid separation 
corresponds to AQ = 1. The torque-free boundary condition, 
if assumed to apply for r < r m j n , implies Ig/rg = ^1/2/ r i /2 • 
All other fluid variables X were simply copied into the guard 
cells, Xg =X\/2, and were subsequently rendered thermody- 
namically consistent. This simple prescription approximates 
free inflow (toward smaller r) across r m ; n . The guard cell val- 
ues for other fluid variables are assigned ignoring curvature of 
the coordinate mesh and formally violate conservation laws at 

The mass of the black hole Mbh was initialized with the 
mass of the initial stellar model contained within r mm . The 
black hole mass was evolved by integrating the mass crossing 
the boundary at r = r^n, 



dM, 



BH 



dt 



= {-Airr pv r ) r=rm 



(49) 



The sum of the mass of the black hole and the mass contained 
on the computational grid remains constant to a high level of 
precision throughout each simulation. 

2.8. Assessment and Tests of the Code 

We conducted tests of internal energy conservation, angu- 
lar momentum transport, and spatial resolution convergence. 



The time-integrated equation for the conservation of internal 
energy in absence of nuclear and thermal energy interconver- 
sion in spherical coordinates reads 



^int,tot(Anax) ^int,tot(^min) — ^^^mm^sl 

P d(r 2 v r ) 



-4tt 



where 



dr 



" £2 vise Qu 



(v r pe+F cmw )dt 



r 2 drdt = 0, (50) 



4nr pedr 



(51) 



and r m j n tes t > r mm is a reference radius defining the inner 
boundary of the spherical annulus in which we test energy 
conservation. We have ignored any flow of energy through 
fm-dx, since stellar material does not reach this radius in the 
course of any simulation. 

In Figure n] we utilize equation (50 1 to test the global con- 
servation ofinternal energy in a run with identical test pa- 
rameters to Run 2, except that we disabled nuclear compo- 
sitional change and used a Newtonian gravitational potential 
without relativistic corrections. In the legend, the apparent 
error AE is defined as the absolute value of the difference 
between the left and right hand side or equation (50i. The 
evaluation of the various terms in equation ( |50| was carried 
out in post-processing from cell-centered data recorded in 
At = 0.01 s intervals, which, in retrospect, is prone to the in- 
troduction of various spatial and temporal discretization arti- 
facts not present in the actual simulation. 

We find that the apparent error is < 1% of the largest 
term in equation d50| when calculated for r m i n test = 50 km 
for the time intervaTOs < t < 70s. The apparent error is 
most significant, AE ~4x 10 50 ergs, prior to and during the 
first few seconds after shock passage. The apparent error 
that accrues after the first few seconds following shock pas- 
sage is less than 10 50 ergs. This can be compared to the to- 
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TABLE 1 

Summary of Simulation Parameters 



10' 



FIG. 2. — The absolute value of" the actual (black, solid) and steady state 
analytic (red, dotted, see equation |52| ) mass flow, M, as a function of radius 
in Run 1 at t = 50 s. The deviation from the analytical value is < 5% at radii 
100 km < r < 1000 km. At this time, r shock ~ 7 X 10 4 km. 



tal binding energy change on the simulation grid, which if 
sufficiently large c an im ply a supernova. We calculate this 
energy in Section |3.5| below and find that it increases by 
~ (1.5-2) x 10 51 ergs following shock formation, and a sig- 
nificant fraction (~5x 10 5() ergs) of the increase is accrued 
later than a few seconds after shock formation, when the 
change in the cumulative apparent error is very small. There- 
fore, it does not seem that the apparent energy conservation 
error at the levels seen in the simulations should significantly 
impact the prospects for explosion. We note that in our calcu- 
lations we explicitly transport specific internal energy rather 
than the total energy, by setting the parameter eintSwitch 
to a very large value. We would like to reiterate that it is likely 
that the apparent error is an artifact of post processing and the 
true energy conservation is better. To demonstrate the latter, 
however, one would have to reconstruct the diagnostic energy 
fluxes using the very same interpolation procedure as is per- 
formed within the PPM in FLASH. We anticipate carrying out 
such a test in an extension of this work. 

In steady state accretion, mass accretion associated with 
viscous angular momentum transport should occur at the rate 



-4-7T 



d_ 
Or 



r vp 



Or 



(52) 



In Figure [2] we show M(r) for Run 2 at t = 50 s along with 
the analytic steady-state estimate of equation ( [52] ). In the ro- 
tationally supported region 100km < r < 1000km, the devia- 
tion from the analytical value is < 5%, which lends credence 
to the accuracy of our angular momentum transport scheme. 

Our spatial resolution was chosen such that we resolve the 
innermost, neutrino-cooled region of the disk over several 
zones. One caveat is that we d o not resolve the sonic radius o f 
the flow, an issue discussed in McKinne y & Gammie| ([2002). 
Because we use a torque-free boundary condition, the calcula- 
tion is not subject to spurious viscous dissipation at the inner 
boundary. However, our boundary condition may still influ- 
ence the fluid flow, and it therefore may be more apt to con- 



Run 


Ar min (km) a 


a b 
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sconv 


p ■ c 
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10.0 


0.1 


0.5 


2.0 


6.0 
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0.5 


0.1 


0.5 
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6.0 


3 


10.0 


0.2 


0.5 


2.0 


6.0 


4 


10.0 


0.025 


0.5 


2.0 


6.0 


5 


10.0 


0.1 


0.5 


5.0 


15.0 


6 


10.0 


0.1 


0.5 


0.5 


6.0 
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0.1 


0.5 


1.0 


6.0 
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10.0 


0.1 


0.25 


2.0 


6.0 


9 


10.0 
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0.5 
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a The minimum resolution element size 

b The dimensionless viscous stress-to-pressu re r atio 

c Rotational profile parameter (see equation |47| ) 

d Convective efficiency parameter (see equation |27| ) 

c Convective compositional mixing efficiency parameter (see equ ation |34| ) 

f This run also had additional angular resolution (see Section |2.7[ 

sider values of M, as opposed to, for example, a or £, when 
comparing our work to other simulations. 

Run 1 and Run 2 contained identical hydrodynamic param- 
eters and differed only in spatial resolution. Run 2 was capa- 
ble of one additional level of resolution refinement over Run 

was set to one- 



1 and the parameter 77 described in Section 2.7 
half the value in Run 1, allowing for significantly higher res- 
olution as a function of radius. Figures|5]|4j|5] and[6]show the 
density, temperature, specific entropy, arid the mean atomic 
weight in Run 1, and also in the higher-resolution Run 2 at 
different times. Substantial agreement is seen between the 
low and high resolution simulations. 

2.9. Limitations of the Method 

The primary limitations of the model of collapsar accre- 
tion that we have presented here include: (1) a very ap- 
proximate one-dimensional treatment of the intrinsically two 
and three dimensional flow structure; (2) limited adequacy of 
the Navier-Stokes viscous stress as a model for the magnetic 
stress arising in the nonlinear development of the magnetoro- 
tational instability; (3) no a priori knowledge of the expected 
efficiency of convective energy transport 4; C onv and of compo- 
sitional transport £conv,mix; (4) treatment of nuclear composi- 
tional transformation through relaxation to NSE rather than 
by integrating the nuclear reaction network that would have 
allowed us to make predictions about the nucleosynthetic out- 
put; (5) neglect of ambient compositional stratification and 
nuclear compositional transformation inside convective cells 
in the calculation of the convective heat flux; (6) the lack of 
modeling of the axial relativistic jet and its enveloping cocoon 
that are thought to be present in LRGB sources; and (7) the 
use of a relatively low mass progenitor star, which may or 
may not be able to yield a black hole and an explosion with 
an energy as high as has been inferred in supernovae asso- 
ciated with LGRBs. Overcoming limitations (1) through (6) 
will require much more computationally expensive multidi- 
mensional hydrodynamic and magnetohydrodynamic simula- 
tions. Limitation (7) can be addressed by applying our current 
method to other, more massive stellar models; here, we spec- 
ulate what collapsars in higher mass progenitors may behave 
like in Section |4]below. 

It would be tempting in view of limitation (5) to try to in- 
corporate the effects of compositional stratification ambient 
to convective cells in the convective energy flux, which is 
normally achieved by multiplying the energy flux in equation 
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FIG. 3. — Density in Run 1 (left) and the higher resolution Run 2 (right) at t = 0s (black, solid), t = 15s (red,dotted), t = 25s (green.dashed), and ( = 50s 
(blue, dash-dotted). In the convective region behind the shock front, some waves form due to an instability, developing at late times, that is likely an artifact of 
including mixing length theory convection as an explicit term in the transport equations. The density jump across the shock front is approximately an order of 
magnitude. 
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FIG. 4. — Temperature in Run 1 (left) and the higher resolution Run 2 (right) at t = 0s (black,solid), t = 15 s (red, dotted), t = 25 s (green, dashed), and t = 50s 
(blue, dash- dotted ). Photodisintegration and neutrino emission cool the innermost disk, while nuclear fusion provides additional heating in the post-shock region 
(see Section[3~4). 



|27) with an factor 



dT\ 



p.p 



dfx 



T ds\ 
CpdrJ 



1 1/2 



(53) 



where fi is the mean nuclear mass, and utilizing the Ledoux 
instead of the Schwarzschild criterion (see, e.g., Bisnovatyi- 



Kogan 2001). This would be meaningful as long as the con 
vective eddy turnover time T conv ~ A conv /v conv were shorter 
than the nuclear time scale r nuc < tnse, so that the convec- 
tive cells can be treated as adiabatic before they mix. How- 
ever at radii where Ledoux convection would differ most from 
Schwarzschild convection, namely where the photodissocia- 
tion into helium nuclei and free nucleons is substantial, the 
convective time scale is much longer than the nuclear time 
scale, T conv T nuc . The internal composition of a convec- 
tive cell evolves as it rises, and the associated entropy change 
is a much stronger effect than the variation of the ambient 
composition treated in Ledoux convection. Magnetization of 
the medium may play a role in this regime but its effects are 
poorly understood. Cognizant of the ongoing research on the 
interplay of convection and nuclear burning (see, e.g., |Arnett| 



|& Meakin|201 l| l, we adopt the Schwarzschild model and con- 
sider it but a parametrization of a complex, still to be explored 
physics. 

3. RESULTS 

Nine simulations were carried out to explore sensitivity to 
the resolution of the simulation Ar m ; n , the viscous stress-to- 
pressure ratio a, the stellar rotation ^, the efficiency of con- 
vective energy transport £ conv , and the efficiency of convective 
mixing £ C onv.mix- The values of these parameters in each of the 
simulations are summarized in Table [Tj Among these, Run 
1 can be considered the fiducial modeL Each simulation was 
run for 100 s, except for Runs 4 and 5, where strong numerical 
instabilities associated with our convection scheme prevented 
us from simulating for more than 40 s and 50 s, r espe ctively. 
In what follows, we present the results. In Section [3~T] we ad- 
dress the evolution of the rate wit h wh ich mass accretes onto 
the central black hole. In Section [372] we discuss the nature 
of radial force balance in the fraction of the stellar material 
that has been traversed by the outgoing shock wave, but has 
not accreted onto the black hole and also discuss the mass and 
angular momentum transport in the system. In Section |3.3| 
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FIG. 5. — The smoothed entropy Jsmooth per baryon in Run I (left) and the higher resolution Run 2 (right) at ; = Os (black, solid), t = 15 s (red, dotted), t — 25 s 
(green, dashed), and t = 50 s (blue, dash-dotted) in units of the Boltzmann constant (k%). After fluid comes into radial force balance, a strong entropy inversion 
is observed, giving rise to convection. 




we address energy transport. In Section 3.4 we address the 



FIG. 6. — Mass-weighted average of the atomic mass A in Run 1 (left) and Run 2 (right) at t = 1 s (black,solid), t = 15 s (red,dotted), t = 25 s (green,dashed), 
and t - 50 s (blue, dot-dashed) . Note that at t = 1 s, the iron core has already accreted onto the central point mass. At late times, photodisintegration in the hottest 
inner regions behind the shock front reduces the value of A. Convective mixing is able to dredge up lighter elements; our scheme for nuclear compositional 
transformation does not correctly model the subsequent recombination and freezeout well outside NSE. 

of the runs, we identify the time when the shock first reaches 
radius 10r m ; n = 250 km as the shock formation time f s hock and 
list the shock formation times in column 2 of Table|2] We also 
provide the mass of the black hole at this point, M BH (f s hock). in 
column 8. The black hole mass at the time of shock formation 
was M B H(f s hock) ~ (5.2-5. 5)M Q in Runs l-7and9. In Run 8, 
which was initiated with reduced initial angular momentum, 
the accretion shock appeared later and the black hole mass is 
correspondingly larger. 

After the formation of the shock, the fluid nearest the in- 
ner boundary is rotationally supported and accretes as a re- 
sult of angular momentum transport driven by the viscous 
shear stress. Subsequent to shock formation, the accretion 
rate declines rapidly either promptly or following a short de- 
lay. The typical rapid drop of the accretion rate is by a fac- 
tor ~ 5- 10 (Runs 1, 2, 3, 5, 7, and 9), and this is followed 
by a continued power-law-like decline. By the end of each 
simulation at 100 s, the accretion rate has typically declined 
to - (10" 3 - lO" 4 )M s- 1 (Runs 1, 2, 3, 4, 7, 8, and 9), or a 
factor of 100- 1000 of the pre-shock value. Final black hole 
masses were ~ (6 — 7)M in the simulations with = 0.5, 



nuclear composition of the flow and discuss the limitations 
inherent in our simplified treat ment of nuclear compositional 
transformation. In Section [3~5] we discuss the global energet- 
ics and check whether sufficient energy may be transported 
into a portion of the stellar envelope to produce a supernova. 

3.1. Central Accretion Rate and Black Hole Mass 

Each simulation exhibits unshocked radial accretion of the 
inner, low-angular momentum mass shells of the progenitor 
star through the inner boundary lasting ~ (20 - 30) s at rela- 
tively steady accretion rates of ~ (0.1 -O.2)M s -1 . The cen- 
tral mass accretion rate, black hole mass, and total mass on 
the computational grid as a function of time in each simula- 
tion are shown in Figure [7] The abrupt drop of the central ac- 
cretion rate at ~ (20-30)s is associated with the appearance 
of an accretion shock precipitated by the arrival of the mass 
shells with specific angular momentum sufficient to lead to 
circularization around the black hole. 

In Figure[8]we show the location of the shock r s h oc k and its 
velocity v s h oc k = dr^ oc ^/dt as a function of time. For each 
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FIG. 7. — Mass of the stellar envelope (top, black, solid), mass of the central object (top, red, dashed), mass of the disk (middle) as defined in Section [TT] and 
mass accretion rate through the inner boundary (bottom) in each of the runs. Most of the mass is accreted onto the central object in the first ~ 20s in most runs. 
The disk mass makes up only a small portion of the total remaining, while the rest of the mass exists in a pressure supported atmosphere that may continue to 
feed the disk, or may be potentially unbound by the accretion shock. Plots of the disk mass begin when t = ? s hock- The quick drop in accretion rate seen in most 
of the simulations occurs around the time of shock formation. 
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FIG. 8. — Shock location (top) and velocity (bottom) in each of the runs. The red dashed line shows rdisk, the outermost radius where the acceleration due to 
the centrifugal force is at least 50% the acceleration due to the pressure gradient. In Run 4 and Run 6, the shock stalls and is reinvigorated. Shock velocities were 
typically 2000— 4000 km s _1 . The small fluctuations in the shock velocity are numerical artifacts of the discreteness in our shock detection algorithm. 
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TABLE 2 

Summary of Key Measurements 



Run 


^shock 


M B Hfthock) b 


^unbound 


^bind 




M Fc f 




1 


20.3 


5.4 


6.0 


0.40 


0.31 


0.06 




28 


19.1 


5.2 


6.4 


0.44 


0.36 


0.04 




3 


19.2 


5.2 


5.7 


0.34 


0.28 


0.06 




4 


19.8 


5.3 


4.4 


0.62 


0.29 


0.07 


cn 


5 


20.6 


5.5 


3.1 


0.40 


0.18 


0.04 


*£ 


6 


20.3 


5.4 


0.0 


-0.43 


0.16 


0.03 




7 


19.2 


5.2 


4.4 


0.08 


0.17 


0.09 




8 


34.0 


7.9 


5.9 


0.54 


0.29 


0.02 




9 


19.2 


5.2 


6.8 


0.46 


0.37 


0.03 





io ,a r 



10 18 



a Time at which shock reaches r = 250 km (s) 

b Black hole mass at when the shock reaches r = 250 kms (Mq) 

c Unbound mass at the end of the simulation (Mq ) 

d Total energy in the s tellar material at the end of the simulation 

(10 51 ergss -1 ; see Section [3~5l and Figure flii) 

e To tal kinetic energy of outbound material (10 51 ergss -1 ; see Section 

1 Total mass of newly synth esize d Fe-group elements at the end of the 
simulation (Mq ; see Section [3~4) 

s This run also had additional angular resolution (see Section [2~7) 

and ~ 10M Q in Run 8 with reduced initial angular momen- 
tum £ e = 0.25. 

In Run 4, with a low value of the viscosity parameter a = 
0.025, the shock first made a very slow progress from 300km 
to 2000km during the first 10 s from its appearance. Then, at 
30 s, the shock suddenly accelerated to v s h oc k ~ 5000 km s~ . 
The near-stagnation of the shock can be understood by notic- 
ing that during the 10 s, the neutrino cooling rate matches the 
viscous heating rate; the rapid cooling prevents the central 
entr opy rise and convection seen i n all other runs (see Sec- 
tion |33] below). In Section 4.2 of |Lindner et~ai~| ( |20T0l >, we 
discussed the scenario in which the shock stagnation brought 
about by efficient neutrino cooling prolongs the LGRB central 
engine activity resulting in a longer prompt emission. 

In Run 6, which had convective efficiency £ conv = 0.5, the 
shock stalled at the radius ~ 10 4 km for ~ 5 s before pro- 
ceeding outward. Note that the reinvigoration of the shock is 
solely driven by the convective energy transport, as we do not 
simulate the negligible neutrino energy and momentum depo- 
sition. The stalling and restarting of the shock was reflected 
in a strong variability of the central accretion rate. 

3.2. The Shocked Envelope and Angular Momentum 

Shock passage leaves a shock- and convection-heated, pres- 
sure supported envelope which contains much more mass than 
the disk, consistent with what we saw in Lindner et al. (2010). 
Figure [3] shows that the density in the envelope is an approx- 
imate power-law of radius p oc r" 9 . Figure indicates that 
the temperature is also a power-law T oc r~°. The pressure 
(not shown) is an approximate power law P oc r" 1 8 . The pro- 
files extend inward into the regime in which rotational sup- 
port dominates pressure support. The mass of the rotation- 
ally supported material in the grid, where a cent > I Igseif+gB-nJ, 
promptly following disk formation was typically < 5% of the 
total mass on the grid. Most of the mass on the grid was in 
the pressure supported atmosphere seamlessly connecting to 
the disk. The mass of the disk in each simulation is shown 
in Figure [7] In some of the runs, certain variability is seen 
in the disk mass over the first few seconds of disk formation. 
Afterwards, the disk mass in each simulation declines mono- 
tonically. In most runs (1-4, 7-9) the disk mass declines to 
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FIG. 9. — Specific angular momentum, t, as a function of radius in Run 1 
at t = 0s (black, solid), t = 15s (red, dotted), t = 25s (green, dashed), and 
/ = 50 s (blue, dot-dashed). The initial rotational profile for the stellar model 
shows large spikes at compositional boundaries. Early in the simulation, low 
angular momentum material is quickly accreted. 



Mdisk ^ 10~ 5 Mq by the end of the simulation, while in Run 6, 
the mass at the end of the simulation is somewhat larger but 
still very small, M disk ~3x IO^Mq. 

Specific angular momentum as a function of radius is 
shown in Figure [9] In the initial angular momentum profile 
of the model, compositional boundaries coincide with discon- 
tinuities in the profile, but in 16TI these occur only at mass co- 
ordinates that are accreted directly onto the black hole, prior 
to the initial circularization. The angular momentum pro- 
file of the mass shells remaining at initial circularization is 
monotonically increasing and most of the remaining mass has 
nearly the same angular momentum, ~ (1— 2) X 10 17 cm 2 s" 1 Q 
This implies that the shocked atmosphere has nearly uniform 
specific angular momentum everywhere except at the radii 
where the time scale on which the viscous torque transport an- 
gular momentum is shorter than the time since circularization. 
At (25-50)s, there is a mild, sub-Keplerian inward downturn 
in £(r) at r < 1000km. Angular momentum transport is too 
slow within the initial ~ 100 s to affect the radii > 10 4 km. 

The mass accretion rate as a function of radius in Run 1 at 
t = 50 s is shown in Figure|2] The accretion rate is independent 
of radius for r < 2000 km, which is the radii where the angular 
momentum profile has relaxed to a viscous quasi-equilibrium. 
The analytic expectation, given in equation ( 52 1, is shown as 
well. Figure [10 shows the radial velocity v r , angular velocity 
and Keplenan velocity as a function of radius at t = 18 s, 
just as material begins to circularize outside of the black hole, 
and at t = 30 s, after an accretion disk has formed. At t = 18 s, 
the velocity reaches the Keplerian value at the innermost 
radii. 

Throughout the simulations, we tracked the value of our 
estimate of the vertical (z-directed) press ure s cale height-to- 
radius ratio Hjr, as described in Section 



2.6 



When the es- 
timated ratio is below one-half, this indicates that in two di- 
mensions, the flow should be disk-like, and when H z /r <C 0.5, 

2 Stellar models exist in which nonmonotonicity is pronounced. T his can 

produce an interesting variability of the centra l accretion rate (e.g., |L6pez-| 

ICamara et al. 2010 Perna & MacFadyen 20T0l. 
I . I II I I 
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FIG. 10. — The absolute value of the radial velocity, v r (black, solid), Keplerian velocity including pseudo-relativistic corrections and self gravity (red, dotted), 
and (green, dashed) as a function of radius in Run 1 at t = 18 s (left) and at t = 30 s (right), just as material begins to circularize. Notice that the rotational 
velocity is approaching the Keplerian velocity at the inner radii. Once material has become rotationally supported, there is a dramatic drop in v r . At radii 
4000 km < r < r s hock. the radial velocity is positive, indicating an outflow. 
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FIG. 1 1 . — The value of the corr ection factor H mod applied to the temperature to account for the possibility of a reduced vertical scale height in, from left to 
right, Runs 1, 4, and 6 (see Section [T6) . This correction factor does not drop below ~ 0.77 in any of the simulations, indicating that geometric thinning of the 
accretion flow is a relatively weak effect. The correction is only applied in regions where a ccn t > 0.5a pl - es , which occurs only for r < 1000 km; thus H raod = 1.0 
for r > 1000km. 



the flow is a geometrically thin disk. We found that H z /r 
is below 0.5 but is still always above a minimum of 0.3 ev- 
erywhere, except in Run 4, which had the lowest viscosity. 
The disk-like radii where the vertical pressure scale height- 
to-radius ratio is below one-half are r < 200km immediately 
following circularization and shrink to r < 100 km by the end 
of the simulation. In Run 4 with a reduced viscous stress-to- 
pressure ratio a, neutrino cooling drove the disk to be geo- 
metrically thin, where H,/r < 0.3 in the inner r < 500 km. In 
the innermost zone in Run 4, H z /r = 0.1 at t = 20s, the lowest 
seen in any simulation. By t = 35 s, no thin disk is present. 
In Figure 1 1 we show the value of S mo( j defined in equation 



44 1 throughout the simulation in Runs 1, 4, and 6; it does 
not drop below ~ 0.77. Only in Run 4 is a genuinely thin 
accretion disk present, and there it is limited to small radii. 



The outer radius of the thin disk decreases as the neutrino lu- 
minosity drops (see Section 3.3 1. We attribute the observed 
moderate thinning of the accretion flow to the cooling of the 
flow by the photodisintegration of helium nuclei into free nu- 
cleons, and in Run 4, the additional contribution of neutrino 
cooling is also significant. 

3.3. Energy Transport 

To understand the energetics of the accretion flow in a col- 
lapsar, we need to consider the transport of mechanical, ther- 
mal, and nuclear binding energy, as well as the loss to neu- 
trino emission. Before turning to energy transport, we discuss 
the neutrino losses, which turn out to be not significant in the 
regime we consider. 

The integrated neutrino luminosity is dominated by the 
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FIG. 12. — Neutrino luminosity integrated over the e ntire computational 
domain in representative Runs 1, 4, and 6 (see Section |2~4) . The peak in 
luminosity occurs shortly after the formation of the accretion shock. Note 
that we do not capture neutrino emission from the region r < 25 km, where 
much neutralization and peak neutrino luminosity is expected to occur in the 
first few seconds after the formation and collapse of the iron core. The total 
neutrino luminosities integrated over the entire simulation in Runs 1, 4, and 
6 were 2.7, 419.3, and 83.9 X 10 5 ' ergs, respectively. 

emission from the inner ~ 100 km. The luminosity as a 
function of time in the representative Runs 1, 4, and 6 is 
shown in Figure [12] In simulations with a = 0.1, neu- 
trino luminosities integrated over the entire computational 
domain peaked immediately following shock formation at 
~ (1 -200) x 10 51 ergss" 1 . The peak luminosity lasted any- 
where from less than a second in the runs with high peak 
luminosities to a few seconds in the runs with low peak lu- 
minosities. After the peak, the luminosity decays first very 
rapidly until it has dropped to ~ 10 50 ergss" 1 , and then con- 
tinues to decay approximately exponentially by several orders 
of magnitude to settle at ~ (10~ 6 - 10~ 5 )ergss _1 after ~ 50s. 
The sharp luminosity peak is an artifact of the abrupt nature 
of shock formation in our 1.5D dimensional treatment and 
is probably not physically significant. The total energy that 
would be deposited by an absorption of the emitted neutrinos, 
which we do not calculate, is negligible. 

Now turning to energy transport, we examine the radial 
transport of all forms of energy, the thermal and kinetic en- 
ergies, the nuclear binding energy, and the gravitational po- 
tential energy. The gravitational potential energy is a non- 
local functional of the mass distribution. However, ignoring 
relativistic effects, one can define the gravitational potential 
energy per unit volume to be /o($bh + |$s)> where $ BH is the 
gravitational potential of the black hole which we define via 
$BH(r) = X°° gBH(r')dr' with g B H given in equation ^ and 
$ se if is that of the self-gravity of the star. Then, pv, $, where 
® = ®bh + $seif, can be interpreted as the flux of gravitational 
energy advected by the fluid, but one must additionally in- 
clude the flux of gravitati onal energy tran sported by self grav- 
ity (see, e.g., |Binney & Tremaine|2 008 Appendix F), which 
equals 



grav.self — 
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Or 



Or 



(54) 



FIG. 13. — The smoothed entropy s smo oth P er baryon in units of the Boltz- 
mann constant (&b) at various times in the low-viscosity Run 4 (cf. Fig. [3). 
Unlike in the other runs, here a negative specific entropy gradient is not seen 
immediately after the fluid comes into radial force balance. Even by t = 30 s 
the fluid is still stable against convection; neutrino cooling prevents the early 
rise of convective instability. Once the neutrino luminosity begins to drop 
around t = 33 s, entropy in the post-shock region begins to rise, bringing about 
strong entropy inversion. By the end of the simulation, Run 4 has the largest 
value of entropy seen in any of the simulations. 

This term is significant only in the outer envelope of the star. 
The rate with which the sum total of these energies is trans- 
ported radially is given by 



E=4nr z 



pv, 



P 1 , 1 I 2 

e+£nuc+-+xV,+ -^+4> 

p 2 2r 



+ F Brz 



grav,self 



on 



(55) 



where the convention is such that E > implies the transport 
of positive energy outward, opposite from the convection em- 
ployed in the definition of the mass accretion rate M. Here, 
-pv£d£l I dr is flux of energy transported by the viscous stress. 

Specific entropy as a function of radius at several times in 
Runs 1 and 2 is shown in Figure E\ After the formation of the 
accretion shock and the rotationally supported flow, viscous 
dissipation heats the fluid, thus producing a negative entropy 
gradient in the shock downstream. The negative specific en- 
tropy gradient extends almost to the shock front, and thus the 
energy injected at small radii can travel to raise the entropy 
of the entire post-shock region. Figure 13 shows the specific 
entropy in Run 4. Here, the high neutrino luminosity after 
the accretion shock has formed keeps the entropy in the post 
shock region relatively low. For the first ~ 10 s after shock for- 
mation, no specific entropy inversion is seen, and the fluid is 
stable against convection. When the neutrino luminosity be- 
gins to drop around t « 33 s, the entropy rises, a negative spe- 
cific entropy gradient appears in the post-shock region, and 
convection starts transporting the viscously dissipated energy 
outward. i 

■ plots the net transport rate E and the various con- 
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stituent terms in Run 2 at t = 30 s; the radii and other observ- 
ables quoted in the remainder of this section will be specific to 
this particular simulation snapshot and will vary across differ- 
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FIG. 14. — The total and partial energy transport rates in Run 2 at f = 30 s (see Section |3.3| and equation [55]). The curves show E (black), the enthalpy advection 
rate 4nr 2 v r (p£+P) (red), the kinetic energy advection rate 2irr 2 v r p(v 2 + £ 2 /r 2 ) (green), the gravitational potential energy transport rate 4nr 2 (pv r Q + F grav jSe i t ) 
(dark blue), the rate of energy transport by the viscous stress -Airr 2 pvldQ. / dr (pink), the rate of thermal energy transport by convection 4irr 2 F C om (green), the 
nuclear energy transport rate associated with convective compositional mixing 4nr 2 F nuc m ; x (gray), and the nuclear binding energy advection rate 4irr 2 v r p£ mc 
(orange). Negative values are indicated by dotted lines. 



ent simulations and different times within a simulation. Ap- 
proximate radial independence of the energy transport rate, 
dE/dr « for 200km < r < 4000 km, where the transport 
rate is positive E w 10 50 ergss _1 > 0, is indicative of quasi- 
steady-state accretion. At larger radii, r > 5000 km, where the 
inner inflow gives way to an outer outflow — a precursor of the 
brewing explosion — no quasi-steady state is present and the 
fluid variables evolve on the dynamical time in the wake of 
the expanding shock. At small radii, r < 100 km, where one 
expects a steady state, the curve E(r) exhibits a small posi- 
tive gradient, as well as a sawtooth consistent with that seen 
in the accretion rate M(r). The constancy of the plotted en- 
ergy transport rate is contingent on an accurate cancellation 
of the other transport terms. We suspect that the observed 
nonconstancy is arising from relatively small inconsistencies 
in the discretization or gravitational source terms in FLASH 
and in the calculation of the gravitational energy during post- 
processing. 

At r < 1000 km, the inward advection of thermal and ki- 
netic energy dominates over the outward transport by convec- 
tion and the viscous stress. Therefore, the inner most flow is 
an advection-dominated accretion flow (ADAF; Narayan & Yi| 



1994| [7995] jBlandford & Begelman|1999| . At r > 1000km, 
the outward transport of thermal energy by convection domi- 
nates the inward transport by advection and this region is thus 
an convection-dominated accretion flow (CDAF, see, e.g., 
|Igumenshch ev et "aL]|2000| [Blandford & 
Inward nuclear binding energy advection 



Stone et al. 1999 



Begelman1 f2604p 
and nuclear compositional mixing both act to transport the to- 
tal energy outward if one counts the negative nuclear binding 
energy in the total energy budget. Convection transports en- 
ergy from the ADAF-CDAF transition radius where the mag- 
nitude of the enthalpy advection flux ~ \v r (pe + P)\ equals 
the convection flux F com to the shock radius r s h oc k- In Fig- 
ure [14] the former is at t"adaf ~ 1000km and the latter is at 
'"shock ~ 3.5 x 10 4 km. In Run 1 and similar runs, tadaf in- 
creases very slowly from ~ 1000 km to ~ 2000km from shock 
formation until t = 100s. In Run 4, the radius is ~ 200km 
throughout the simulation, and in Run 6, the radius grows 
from ~ 5000km to over ~ 10 4 km in the course of the sim- 
ulation. 

We suspect that the location of tadaf determines the 
amount of energy that can be carried to the shock front, and 
that in turn, the ADAF-CDAF transition radius is primarily 




FIG. 15. — The abundances of the common elements summed over all isotopic species at (left to right) t = 0, 15, 25,50s in Run 1. 



a function of the convective efficiency £ conv and the viscous 
stress-to-pressure ratio a. Simulations with larger values of 
Cconv resulted in stronger shocks and larger amounts of un- 
bound material. The convective compositional mixing param- 
eter £ C onv.mix had little effect on the final outcome of our simu- 
lations. Turning to the viscosity parameter, Run 3 with a large 
values of a produced a somewhat less energetic shock with 
less unbound material at the end of the simulation. Run 4, 
the simulation with the lowest viscosity, produced the most 
energetic explosion, even in the presence of more pervasive 
cooling by neutrino emission and by photodisintegration in 
the 1ow-q regime; these flows are denser and hotter (s ee, e.g., 
I Popham et al |1999[|Cheh & Beloborodov|2007[ l. In] Milosayl-| 
[jevic et al.| ( |2012| ) we show that ^adaf is expected to be smaller 
under low viscosity conditions because the advective luminos- 
ity is proportional to v, , which is proportional to a. This trend 
is reproduced in the present simulations. 

3.4. Nuclear Composition of the Flow 

Our simplified treatment of nuclear compositional transfor- 
mation, which entails relaxation to NSE on a temperature- 
dependent time scale, is designed to track the impact of nu- 
clear photodisintegration and recombination on the thermo- 
dynamics of the flow. However, it does not allow the compu- 
tation of the ultimate nucleosynthetic output in the presence of 
out-of-NSE burning. Thus, the results presented here can only 
be understood in the light of the limitations of the method j^| It 
is also worth recalling that we do not calculate neutronization 
that could modify the proton-to-nucleon ratio Y e . 

In the hottest, innermost accretion flow, photodisintegration 
of heavy nuclei saps energy that could otherwise be trans- 
ported by convection to larger radii to energize the shocked 
envelope. However, once the nuclei are broken down, convec- 
tive mixing can dredge up free nucleons to larger and cooler 
radii, where they can recombine and heat the fluid locally. 
Figure|6]shows the mean atomic mass A as a function of radius 
at various times in Run 1 and Run 2. The mean atomic mass 
drops below_4 in the innermost (200-300)km. The positive 
gradient in A seen in portions of the convective region would 
in the Ledoux picture enhance the convective energy flux, but 
our Schwarzchild treatment o f con vection does not capture 
this effect. We argue in Section|Z9]that since the nuclear time 



Metzeer ](2012} modeled accretion disks associated with the mergers of 
white dwarfs and neutron stars or black holes. He found that nuclear pro- 
cesses taking place at temperatures T < 4 X 10 9 K may lead to significant 
heating in the resulting outflows. The relaxation to NSE we employ underes- 
timates the heating due to out-of-NSE nuclear recombination. 



scale is shorter than the convective time scale at radii where 
Ledoux convection implies an enhanced energy transport, the 
nuclear compositional transformation inside convective cells, 
not consider in the Ledoux treatment, should dominate. Lack- 
ing a theory of convection in this regime we adhere to the 
simpler Schwarzschild parametrization. 
In Figure 15 we show the mass-weighted abundances of 



the most common elements in our simulations in Run 1 at 
various times. By t = 30s, again, the inner 200km is made up 
almost entirely of free nucleons in nearly equal portions, as 
Y e 0.5 everywhere. The effect of convective compositional 
transport of the reprocessed nuclear species — the free nucle- 
ons, helium, and iron — from the hot innermost accretion flow 
is seen in the power law tails extending to near the location of 
the shock in the right panels. 

Although in our calculations we do not allow the evolu- 
tion of Y e , we can still specula te about the effects of neutron- 
ization. |Chen & Beloborodov| ( |2007] l computed the structure 
of time-independent accretion disks around Kerr black holes 
including the effects of pair capture and neutronization. In 
their models with a = 0.1, the same as our fiducial viscous- 
stress-to-pressure ratio, at the radii where p ~ 10 7 g cm" 3 , cor- 
responding to the density in the innermost disk in our simula- 
tions, they find Y e sa 0.5, the same as in our non-neutr onizing 
treatment. In their models with a = 0.01, however, |Chen &| 
Beloborodo~v]( |2007| find that at densities corresponding to the 
innermost disk in our simulations, significant neutronization 
was in effect and Y e dropped well below neutron-proton equal- 
ity. It is therefore possible that in the very innermost regions 
of the disks in our Run 4 with a low viscosity a = 0.025, the 
true value of Y e should be lower than we assume. This would 
modify the abundances and thermodynamics of the portion 
of the flow in NSE. The key question of consequence for 
the viability of the mechanism we propose for the produc- 
tion of luminous supemovae is, will the neutron-rich material 
pollute larger radii and drive a tendency toward the synthe- 
sis of iron instead of 56 Ni? Because this neutronization only 
seems to be most significant in the hottest innermost regions, 
where neutrino cooling is efficient and the flow is predomi- 
nantly rotationally supported, it is possible that most of the 
neutron-enhanced material would be advected into the black 
hole. This is a quantitative question that can be answered only 
with multidimensional simulations. 



In Figure 16 we show the location of r N $E> the largest ra- 
dius at whicn nuclear compositional transformation, in the 
form of gradual convergence toward NSE, is taking place in 
our Run 1 . This is the only region in which our calculations 
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FIG. 16. — Location of ^nse, the largest radius at which nuclear compo- 
sitional transformation, in the form of gradual convergence toward NSE, is 
taking place in our Run 1 ; this radius is closely approximated with the radius 
where the temperature crosses ss 3 X 10 9 K. 



will capture nucleosynthesis and photodisintegration. Effec- 
tively, it is the radius at which T > 3 x 10 9 K. Early in the sim- 
ulation, the hot stellar core accretes into the black hole, and 
'"nse quickly drops from « 3000 km to « 800 km. After the 
shock forms, additional heating in the shock and by viscous 
dissipation rapidly increases r N sE and it peaks at s» 4000 km. 
As the density and temperature drop and the viscous heating 
rate decreases, the inner regions of the star begin to cool, and 
tnse declines again. 

In Figure [TTj we show the total mass of various nuclear 
species in the entire simulation domain as a function of time. 
The most notable change is the dip in the mass of iron-group 
elements. The initial decrease in in the iron group mass is the 
accretion of the core of the star onto the black hole. After 
shock formation, a rapid increase in the amount of free nucle- 
ons is seen, in addition to production of additional helium and 
iron group elements. In Table |2j we show the total amount of 
newly fused Fe-group elements present at the end of the sim- 
ulation, which fall in the range of O.O2M Q -O.O9M . Since 
we do not calculate the out-of-NSE burning in convectively 
dredged up material, at least a fraction of the extended he- 
lium tail seen in Figure[l5]could be expected to burn into iron 
and thus the iron group mass in Figure[l7]and Table[2]can be 
interpreted as a lower limit. 

3.5. Prospects for Explosion 

The total thermal and mechanical energy, which we also 
refer to as the total binding energy, present on the grid was 
computed for each simulation via 



^bind — 



and is shown in Figure 18 A positive binding energy indicates 
a potential for explosion. Runs 1, 2, 3, 4, 5, 8, and 9 acquired 
a positive total binding energy £bind ~ (3.5-6.2) x 10 50 ergs 
by the end of the simulation, indicating the potential for ex- 
plosion. These runs have £ conv > 2 in common. Run 7 
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FIG. 17. — Mass of the most common elements summed over isotopic 
species in Run 1 as a function of time. Here, Fe represent all the isotopes 
of the iron group. At the start of the simulation, there is a large dip in the 
amount of iron group elements, due to the accretion of the iron core. After 
shock formation, ~ 0.05 Mq of iron group elements are synthesized. 

with £ conv = 1 reached marginally unbound condition with 
£knd ~ (0.5-1) x 10 50 ergs. Run 6 with £ conv = 0.5 remained 
gravitationally bound throughout the entire simulation. 

The total unbound mass in each simulation was calculated 
by summing the masses of any fluid element with a positive 
value of £bind an d a positive radial velocity. The total un- 
bound masses in each simulation defined by this diagnostic 
are shown in Table [2] This criterion does not take into ac- 
count any interaction that unbound and bound material might 
have subsequent to the measurement. This criterion also ne- 
glects future energy gain or loss from nuclear processes. The 
location and velocity of the outward moving accretion shock 
is shown in Figure [8] In exploding models, the typical shock 
velocities were (2000 -4000) km s" 1 . In Run 4 and Run 6, 
the shock stalled or slowed, and later restarted once or sev- 
eral times. In Figure [19] we show the evolution of various 
Lagrangian mass coordinates in Run 1 throughout the simula- 
tion. Once the shock moves beyond ~ 4000 km, the infalling 
material obtains a positive velocity once it reaches r s h oc k- 

These results suggest that a high convective efficiency is 
required for sufficient transfer of energy from the inner accre- 
tion flow to the envelope to unbind the envelope. Simulations 
with higher values of £ conv had relatively larger shock veloc- 
ities and amounts of unbound mass at the ends of their sim- 
ulations, and Run 4 with a lower a had the largest value of 
£bind- Run 8 with reduced initial specific angular momentum 
produced an explosion comparable to that seen in the fiducial 
model. 

4. OBSERVATIONAL SIGNATURES AND PROGENITOR TYPES 

Modeling of light curves and spectra of supernovae asso- 
ciated with LGRBs has yielded information about the nature 
of these explosions. The mass of 56 Ni present in the super- 
nova ejecta is easily estimated from the light curve by fitting 
simple radioactive decay models. The velocity of the ejecta 
is inferred from observed line widths. Then, in a standard ap- 
proach, the kinetic energy and mass of the ejecta are derived 
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FIG. 18. — The total binding energy on the simulation grid in Runs 1, 
and 2—9. Simulations in which the convective efficiency factor, £ C onv > 0.5 
reached a positive total binding energy by the end of the simulation. The 
large, brief dips in energy visible in some simulations around the time of 
shock formation are primarily due to cooling by neutrino emission. The large 
dip seen in Run 4 is also predominantly due to neutrino cooling and has a 
minimum of iJbind = —1.5 X 10 52 ergs. 

by comparing the inferred 56 Ni yield with that implied by one- 
dimensional hydrodynamic models in which a spherical shock 
wave is introduced by hand (by the action of a hypothetical 
"piston") into the progenitor's envelope]^] This approach re- 
lies on the hypothesis that 56 Ni is produced at the shock front, 
for which the shock must be very strong. 

Our simplified simulations suggest an alternative scenario 
in which nucleosynthesis takes place in the accretion flow in 
the interior of the star, similar to the wind-nucleosynthesis 
models (e.g., |Beloborodov|2003|rPruet et al.|2003| [20041 |Na^1 



In "piston" and "thermal bomb" models, an explosion is mimicked by 



in jecting a large kinetic or thermal energy into a narrow shell over a relatively 
short < 1 s time interval (e.g., Jones et al. 1981 Woosley & Weaver 1986 
Arnett|198 /||S higeyama et al.|198 /|| Woosley et al.|1988||lhielemann et all] 
lWOHWo osley & Weaver|lW5||Llmongi &' Chiefh 2003 , Chiefh & Limongi 
2004[|Young & ]Vryer|2^07||l''ry"r et al . 2008 Kasen & Woosley 200 9||Maeda| 
& lbminaga|2U09[|Joggerst et aI.|20lU||L>essart et al.|201 If. 
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FIG. 19. — The evolution Lagrangian mass coordinates in Run 1 (black, 
solid). The location of r s hock is also shown (red, dashed). The bifurcation 
between an inner inflow and an outer outflow occurs at r ~ 4000 km. This 
corresponds to a mass coordinate of ~ 5.6Mq. 



gataki ^et^[2q06l|Surman et al.|[2006l |Maeda & Tom inaga 



2009| |Metzger 2012) . Our results, however, suggest that the 
accretion flow is long lived, lasting tens or hundreds of sec- 
onds or longer, and so the nucleosynthesis can be sustained at 
lower densities than in the wind models, and its products can 
be delivered to the envelope by vigorous convection. 

We also find that a supernova-like shock wave may be pow- 
ered by the sustained input of accretion energy, without ener- 
gization by neutrino energy deposition. The dynamics of the 
accretion-energy-powered shock wave is fundamentally dif- 
ferent from that powered by a piston. It remains to be ex- 
plored whether the accretion scenario will call for a modifi- 
cation of the standard approach to modeling the light curves 
and spectra of the supernovae that could be yielding black 
holes. We are thus somewhat reluctant to compare our re- 
sults directly with observational inferences obtained with ex- 
isting supernova models. Previous work has attributed kinetic 
energies of ~ (2 - 50) x 10 51 ergs to supernovae associated 
with LGRBs (see, e.g., |Woosley & Bloom|[2006{ |Hjorth &| 
Bloo m|2011| and references therein). Our models come short 
of these energies, but they are consistent with the low energy 
end among the more typical T ype lb and Ic supernovae (see, 
e.g., Table 4 in Hamuy 2003). Unfortunately the simplified 
treatment of nuclear compositional transformation does not 
allow us to predict the 56 Ni synthesized in our models. We 
can only say that supernovae powered by collapsar accretion 
should exhibit a high degree of mixing of hydrostatic and ex- 
plosive elements. 

The shock velocities at t = 100 s, when the shock is still in 
the interior of the envelope, in the models that achieve explo- 
sion, are v s h oc k ~ 4000 kms -1 . This is a half or smaller fraction 
of the commonly cited values for shock velocities measured in 
the observed supernovae. Of course, leading to the breakout 
of the stellar surface, the shock accelerates as it traverses the 
negative density gradient. The mass-weighted rms free expan- 
sion velocity inferred from the total energy at the end of the 
simulation is v FE ~ (2£ , bind /M unbound ) 1 ' /2 ~ 3000(£ b i n d/0.5 X 
10 51 ergs) I / 2 (M un bound/5M )" 1/ ' 2 kms" 1 , again lower than usu- 
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ally quoted for the observed supernovae. 
Our initial model of choice was the M stal 



14M Q Wolf- 



Ray et model star 16TI of Woosley & Heger (2006), evolved 
to pre-core-collapse from a 16M Q main sequence progeni- 
tor. The model 16T I is commonly used in LGRB investi- 
gat ions ( e.g., Morso ny et al.|[2007| |2010| |L6 pez-Ca mara et 
aD|2009l 120101 ILazzafi et al. 120101 1201 lalbl ILindner et al. 



2010 Nagakur a et al.||201 la|b[ l, "but it has been suggested 
that the progenitors or supernovae with confirmed associa- 
tion with LGRBs must be associated with the core collapse 
of more massive stars, perhaps wi th masses in the range of 
(25-6O)M or higher (e.g., |Podsiadlowski et al.|2004|[Smartt| 
2009 and references therein). However, predictions regard- 
ing the nature of the final remnant in such explosions are sen- 
sitive to the highly-debated details of the explosion mecha- 
nisms in core collapse supernovae. Observational studies of 
the spatial distribution of Type Ic supernovae and GRBs in 
galaxies suggest that the respectiv e progenitors sh ould be at 
least ~ 25 M and ~ 43 M (Raskin et al.||2008| see, also, 



ILarsson et 



iviq aim ~ h 

al.|2007| >. It 



is of interest to note that simulations 
of neutron-star-powered explosions have been successful only 
in the lowest mass progenitors. The accretion powered mech- 
anism we propose will operate in more massive progenitors 
that produce black holes. It is reasonable to speculate that in 
progenitors more massive than in our model, or with differ- 
ent internal structure, the explosion energies would be much 
higher than we find, more in line with the high energies of 
the LGRB supernovae. The long term accretion in massive 
collapsar progenitors deserves further study. 



5. CONCLUSIONS 

We have conducted a series of hydrodynamic simulations of 
the viscous post-core-collapse accretion of a rap i dly ro tating 
— 14M Wolf-Rayet star of [Woosley & Heg er (2006} onto 
the central black hole that we assumed was in place at the 
beginning of the simulation. The spherically-symmetric sim- 
ulations with rotation were carried out for 100 s and resolved 
the radii down to 25 km, including the range of radii where the 
collapsing stellar material circularizes around the black hole. 
The simulations tracked cooling by neutrino emission and the 
relaxation to nuclear statistical equilibrium in the hot inner 
accretion flow. The simulations also tracked convection and 
convective compositional mixing in the mixing length theory 
approximation. Finally, the simulations tracked viscous an- 
gular momentum transport and the associated heating in the 
flow. To explore the sensitivity to model parameters, we var- 
ied the initial angular momentum profile, convective energy 
transport and compositional mixing efficiencies, and the vis- 
cous stress-to-pressure ratio a. Our main findings are as fol- 
lows. 

Lacking sufficient angular momentum to be rotationally 
supported around the black hole, the inner mass coordinates 
of the stellar envelope accrete unshocked onto the black hole. 
At t ~ 20 s, or later with reduced initial angular momentum, 
the first mass shell able to circularize around the black hole ar- 
rives. Once material becomes circularized, an accretion shock 



forms as the mass shells in near free fall interact with the ro- 
tationally supported material. 

This shock front travels outwards, leaving behind a mostly 
pressure supported, shock heated, convective stellar enve- 
lope. Only a very small fraction of the mass is predominantly 
rotationally supported; the rotationally supported, geometri- 
cally thick disk connects smoothly to the pressure supported, 
shock-heated envelope. The structure and energetics of the 
flow are governed by accretion mechanics. The energy dis- 
sipated by the viscous stress at the innermost radii, the radii 
smaller than some critical r AD AF, is advected into the black 
hole. The innermost flow is thus an ADAF. At larger radii, 
convection transports the dissipated energy outward, into the 
stellar envelope and toward the expanding shock front, imply- 
ing that the outer flow is a CDAF. The delivery of energy from 
~ tadaf to the envelope proceeds for tens of seconds, and the 
total energies delivered are sufficient to produce supernovae, 
albeit not as energetic as the ones inferred in association with 
LRGBs. 

We found that the final energy deposited into the enve- 
lope strongly depended on the efficiency of convective energy 
transport and depended somewhat on the viscous stress-to- 
pressure ratio a. These two parameters strongly influence 
the location of the ADAF/CDAF transition, as we have ex- 



plored with crude analytical arguments in Milosavljevic et al. 
(|2012|). The low-a model has a hotter disk with more per- 



vasive cooling by neutrino emission and nuclear photodis- 
integration, sapping energy from the final explosion. How- 
ever, the low-a disk also has an ADAF/CDAF transition at 
a smaller radius, potentially allowing for a higher convective 
luminosity. 

For sufficiently high convective efficiencies, the stellar en- 
velope was capable of obtaining positive total thermal and 
mechanical energies ~ 0.5 x 10 51 ergs, shock velocities ~ 
4000km s -1 , and unbound masses ~ 6M . We suggest that 
the accretion powered mechanism, which is distinct from and 
possibly mutually exclusive with the standard neutron-star- 
powered "delayed-neutrino" mechanism, could explain low 
luminosity Type lb and Ic supernovae, but multidimensional 
study is needed to pin down the true efficiency of convective 
energy transport and to estimate the expected 56 Ni yield. 
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